The fate of the resonating valence bond in graphene 
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We apply a variational wave lunction capable of describing qualitatively and quantitatively the so 
called "resonating valence bond" in realistic materials, by improving standard ab initio calculations 
by means of quantum Monte Carlo methods. In this framework we clearly identify the Kekule and 
Dewar contributions to the chemical bond of the benzene molecule, and we establish the correspond- 
ing resonating valence bond energy of these well known structures (~ O.OleV/atom). We apply this 
method to unveil the nature of the chemical bond in undoped graphene and show that this picture 
remains only within a small "resonance length" of few atomic units. 



Since the recent experimental isolation of 2D graphene 
layers [lj, there has been a renovated interest in the elec- 
tronic properties of graphene. On the other hand the res- 
onating valence bond (RVB) theory was proposed several 
years ago by Linus Pauling [2] and its successful applica- 
tion to aromatic compounds containing the benzene ring, 
has immediately raised the question whether this fasci- 
nating theory remains meaningful in graphene, which can 
be viewed as a two dimensional realization of Carbon 
rings in a honeycomb lattice. 

Graphene is a subject of intense studies, also because 
its peculiar band structure implies a vanishing density of 
states at the Fermi energy with Dirac cones and non con- 
ventional semimetallic behavior [3| .We also mention that 
the photoemission properties, and a possible opening of 
a gap around the Dirac cones have not been fully un- 
derstood neither experimentally 0] nor theoretically]!, [f| , 
and recently it has been speculated that electron correla- 
tion may play a crucial role in this material [7j, and could 
lead not only to the explanation of this effect but also 
to a rather speculative d + id (room)-high-temperature 
superconductivity upon doping. Generally speaking the 
role of electron correlation in graphene remains highly 
controversial [8], and the attention in the field has been 
renewed by a recent numerical simulation of the Hubbard 
model on the honeycomb lattice|9j- In that work, by us- 
ing an unbiased numerical method, it was shown that the 
ground state of the model could be highly non trivial: an 
insulator, with neither magnetic nor whatsoever broken 
symmetry, namely a RVB spin liquid state. 

In this Letter we clarify the role of RVB correla- 
tions in graphene and other Carbon compounds by us- 
ing a tool[10( for ab-initio calculations based on quan- 
tum Monte Carlo (MC) methods, capable of describing 



rather well the electron correlation in several challeng- 
ing molecules, up to the quantitative description of the 
weak binding in graphite ll| . With this technique we 
can visualize the RVB character of the chemical bond, 
and describe realistically an RVB spin liquid state, with 
the same type of variational wave function that has been 
shown to be rather accurate in model systems, especially 
in spin ones [l2| ■ Since in realistic models that allow 
charge fluctuations, like e.g. the Hubbard model, it is not 
possible to work with a complex wave function without 
breaking time reversal symmetry, we restrict our varia- 
tional freedom to real wave functions, which nevertheless 
allow a very wide class of spin-liquid states. 

We shortly describe the wave function used in this Let- 
ter, as more details have been published elsewhere (see 
e -g- EII and refs. therein). The RVB ansatzfloT] 
|RVB) = J|AGP) (JAGP) is made of a product of a Jas- 
trow factor </, which takes into account the short range 
strong Coulomb repulsion, and the so called antisym- 
metrized geminal power (AGP). A singlet valence bond 
between two electrons of opposite spin is determined by 
a geminal function /. At variance with the usual Slater 
Determinant (SD), where no correlation between oppo- 
site spin electrons is considered, in the AGP all the elec- 
trons are paired with the same geminal. The resulting 
wave function is then antisymmetrized. We parametrize 
/ by using a given number n* of molecular orbitals 
(MOs) as /(r^r 1 -) = YZ n k^k(r^)ipk(ri), where n k 
are variational parameters. The MOs ipk are expanded 
in an atomic basis set and fully optimized by minimiz- 
ing the variational MC (VMC) energy expectation value 
of the full electron-ion Hamiltonian within the Born- 
Oppenheimer approximation [lj]. In all the calculations 
of this work, we have replaced the Is core electrons of 
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the Carbon atom with appropriate pseudopotentials fl6jj . 
which also account for scalar relativistic effects. When 
n* > N/2, with N denoting the number of electrons, 
the wave function has a larger variational freedom with 
respect to the best (lowest in energy) Jastrow SD wave- 
function (JSD)|17|, and is able to improve the descrip- 
tion of the electron correlation, especially when the AGP 
is used in combination with the Jastrow factor. The 
latter is particularly important for the description of 
a spin liquid state and is represented by a weight fac- 
tor J(R) = exp[^ i< ■ u{fi,rj)] over the 3iV-dimensional 



Molecule 


(V)JSD 


(V)JAGP 


(LR)JSD 


(LR)JAGP 


Exp. 


c 2 

C 6 H 6 


5.54(2) 
56.98(1) 


6.33(2) 
57.11(3) 


5.76(2) 
57.11(1) 


6.30(2) 
57.14(1) 


6.30(2) a 
56.62(3) 6 



TABLE I: VMC (V) and LRDMC (LR) atomization energy 
(in eV) of C 2 (AGP primitive basis: 5s5p) and C 6 H 6 (C prim- 
itive basis: 24s22pl0d6,f; H: 3s2p). 



configuration R of the electron positions r^. For the ex- 
plicit form of J, see e.g. Ref. [lj. Provided the two- 
electron function u(fi, fj) decays slowly enough with the 



Molecule 


C 6 H 6 


8C 


16 C 


48 C 


All 


0.118(2) 


0.159(7) 


0.207(4) 


0.18(1) 


7T 


0.101(2) 


0.116(5) 


0.147(8) 


0.15(1) 



distance between the electrons 



|, it is possible to 



TABLE II: VMC contribution (in eV) of all (All) the occupied 
bands and of the n band to the binding energy of C 6 H 6 and 
graphene layers of 8C, 16 C, and 48 C atoms (AGP primitive 
basis: Ils9p7d). 



describe rather well a spin liquid insulator, even when, 
in absence of J(R), the AGP pairing function describes 
a semimetal (for n* = N/2) or a superconductor (for 
n* > n)[l8j]. As discussed in Ref. [3, an appropriate 
choice of n* is crucial to improve the accuracy in the de- 
scription of the chemical bond with the JAGP ansatz: 
n* is the minimum number of MO's that can be used for 
describing a product of independent Hartree-Fock (HF) 
wave functions for isolated atoms. Within this choice 
of n* , both the two-electron functions / and u are ex- 
panded in a basis of localized gaussian atomic orbitals, 
with a method that in principle converges to the com- 
plete basis set limit (CBS), yielding the lowest possible 
energy state compatible with the given ansatz 19] . 

We test our variational ansatz on small Carbon com- 
pounds. We consider the atomization energy of the Car- 
bon dimer and of benzene, computed as the difference 
between the JAGP energy for the entire molecule and 
the JSD energy of the isolated atoms [13] . To compare our 
results with the experimental binding energies we also in- 
clude inner shell correlations and relativistic effects and 
we subtract the zero-point energy. In Table H] we show our 
VMC and lattice-regularized diffusion MC (LRDMC) 
results. The simulations for the Carbon dimer were per- 
formed in Ref. [13 (n*=7). We evaluate inner shell corre- 
lations by comparing the all electron C, energy found in 
Ref. [2l| with the energy found in Ref. HH, where the same 
pseudopotential[16| as for Ref. 1 141 was used. We take spin 
orbit effects from Ref. 0. For benzene (n*=24) we use a 
large basis set close to the CBS-limit within O.OleV/atom 
in all the cases studied. We take inner shell correlation 
and spin orbit effects from Ref. One of the main 
results of our calculation is represented by the sizeable 
energy gain that is obtained by using a large number 
of MOs in the AGP part of the wave function. This 
energy gain is particularly important to get a quantita- 
tive description of the C 2 chemical bond, whereas it is L y ~ m\/3a where a = 1.42 A is the nearest neighbor 



ods/experiments. This is plausible considering that 1) 
the result provided by LRMDC, which is known to im- 
prove the energy estimate [2Cj. is in agreement with our 
best variational ansatz, and that 2) we have never over- 
estimated the well depth by more than O.OleV/atom in 
all the cases studied in the previous work[14]. Anyway, 
our variational ansatz appears to be adequate and en- 
courages us to quantify the RVB energy, which, in the 
present formulation, can be defined as the energy dif- 
ference between the best variational energy found with 
n = N/2 MOs and the one with n = n* > N/2, both 
obtained in presence of J. In Table UH we report the 
contribution of the 7r-band orbitals to the RVB energy 
of benzene and graphene. The tt orbitals yield approxi- 
mately the 80% of the pairing and represent in general 
the most important contribution, as expected. 

To get a deeper insight into our variational calcula- 
tion with the JAGP wave function, we introduce also 
a "valence-projected pairing function" (VPPF), defined 
as fvppF{rt,r±) = J2k>N/2 n k t Pk(rf)tpk(ri)- In the HF 
case of a single SD, fvppF{r^,r^) = 0. Hence, when 
singlet valence bond pairing occurs and is non zero 
even for k > N/2, we can visualize and characterize, in 
real space, the genuine RVB contribution to the chemi- 
cal bond. We can also plot the VPPF restricted to the 
7r-band as a function of r± as done for benzene in Fig. [T] 
Kekule nearest neighbor and Dewar further neighbor cor- 
relations are manifest. Fig. [1] proves the JAGP wave 
function to be a powerful tool for the description of the 
fundamental features of the RVB chemical bond. 

We now discuss the case of undoped graphene. We con- 
sider rectangular supercells L x x L y , with L x = 3na and 



possible that the slight overestimation of the atomiza- 
tion energy for benzene does not depend on the accuracy 
of our total energy estimates, but comes out from the 
previously described corrections taken from other meth- 



Carbon distance and n, m are integers. We use an in- 
creasing number 4nra of C atoms (8, 16, 24, and 48, with 
n, m such that L x / L y ~ 1). These supercells do not sat- 
isfy the 7r/3 rotation symmetry of the infinite lattice 25 1. 
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FIG. 1: (Colors online) Two dimensional plot of the AGP 
pairing function restricted to the molecular orbitals above the 
HOMO. The arrow indicates the reference position r-j- fixed 
on an atom, colored in red for the sake of clarity. 



FIG. 3: (Colors online) Ratio between the s-wave and the 
d-wave weight in the J AGP wave function for the CaCu0 2 
parent compound of the high-temperature cuprate supercon- 
ductor (2x2 supercell). 
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FIG. 2: (Colors online) Two dimensional plot of the AGP 
pairing function for a graphene layer of 48 C atoms restricted 
to the molecular orbitals above the HOMO (VPPF). The ar- 
row indicates the reference atom. 



This helps the system to break rotational symmetries, 
y 2 for a real pairing function, that are 



such as d xy or d 



energetically favored when the expected d + id pairing 
symmetry [26[ characterizes the ground state wave func- 
tion. For each system size, we optimize the JAGP wave 
function, find the VMC energy and the VMC RVB energy 
(by means of correlated-sampling simulations). As re- 
ported in Table HI1 we have also checked the contribution 
of the 7r-band orbitals to the RVB energy gain. In Fig. [2] 
we show the VPPF (restricted to the 7r-band) for the 
largest supercell considered here. Despite the small num- 
ber of atoms, we already see an almost perfect rotational 
symmetry of the VPPF, that is not compatible with d- 
wave pairing. To prove that our method is capable of 



tackling with pairing functions with d— wave symmetry 
we apply our scheme to the CaCu0 2 parent compound 
of cuprate high-temperature superconductors. As shown 
in Fig. [31 in less than 3000 iterations we melt the s-wave 
pairing and are able to detect the correct d-wave symme- 
try of the pairing function. We can conclude, therefore, 
that the RVB chemical bond in graphene is characterized 
by a pairing function with a clear s-wave symmetry. 

Finally, in order to understand the thermodynamic 
properties of graphene, we consider a finite size scal- 
ing of our results. In Fig. 0] we show the energy gain 
due to the s-wave RVB (upper panel) and the ratio 
R = n N / 2+ i/n N / 2 01 the LUMO/HOMO weights as 
a function of the inverse number of C atoms in the su- 
percell. Before discussing this result, we recall what hap- 
pens to the above-mentioned quantities in the absence 
of correlation, i.e. when there is no Jastrow factor in our 
variational ansatz. In such a case, if the ratio R con- 
verges to a finite quantity in the thermodynamic limit, 
the AGP wave function describes an s-wave supercon- 
ductor with true off diagonal long range order. Besides, 
the thermodynamic-limit RVB energy per atom remains 
finite and represents just the condensation energy of the 
s-wave superconductor. In the presence of J instead, a 
different scenario is possible. Indeed, a ratio R > in 
the thermodynamic limit and a finite RVB energy /atom 
denote a spin liquid state with a spin and a charge gap 
in its spectrum. This possibility is compatible with the 
recent Hubbard model results |9|, and may explain also 
the existence of a small gap in the photoemission experi- 
ments, genuine and determined only by the RVB charac- 
ter of the ground state. Due to the computational cost 
increase for larger supercells, it is difficult to obtain an 
accurate thermodynamic limit with our VMC method. 
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FIG. 4: (Colors online) (a): RVB energy per atom for the 
graphene layer, as a function of the number of atoms in the 
supercell (r point), (b): Ratio of the LUMO/HOMO weight 
in the AGP, a measure of the RVB character of the bond. 
Lines are guides to the eye. 



However, clear trends are evident from Fig. |4j In the up- 
per panel, we see that the energy gain of the RVB wave 
function systematically decreases as the system size in- 
creases, apart for the negligible value found for the 24 C 
supercell. The anomaly of the 24 C cluster can be easily 
explained as a shell-effect. Indeed, this cluster should be 
closer to the thermodynamic limit, since it contains the 
so-called K point, the gapless Dirac point in graphene. 
This shell effect does not affect the eigenvalues of the 
pairing function, which instead decrease monotonically 
as the system size increases and reach a very small value 
in the thermodynamic limit (lower panel). If we extrap- 
olate the upper-panel results, omitting the 24 C cluster, 
also the RVB energy per C atom becomes extremely small 
in the thermodynamic limit (smaller than the accuracy 
of the present data) . Both panels thus suggest that the 
semimetal character of graphene should be stable in the 
thermodynamic limit. A small gap could appear in the 
excitation spectrum only if its value was extremely small 
~ O.OleV. We have estimated this value by matching our 
results for the with the ones obtained with an s-wave 
BCS hamiltonian with nearest and next-nearest neighbor 
coupling, describing a Z2 gapped spin liquidQ when cor- 
relation is included by means of an appropriate Jastrow. 

In conclusion we have systematically studied Carbon- 
based compounds from the simplest C 2 molecule to 
graphene layers. We have shown that the RVB character 
of the chemical bond can be depicted in terms of a very 
powerful and accurate wave function that not only im- 
proves the description of the chemical bond but it is also 
capable to show qualitatively new effects induced by the 
electron correlation. We have found clear numerical evi- 
dence that singlet s-wave pairing in graphene be quite ro- 



bust and sizeable up to a small length scale of few atomic 
units. This feature might remain in the thermodynamic 
limit leading to a very small gap in the photoemission 
spectrum or to s-wave superconductivity upon doping, 
effects that can be in principle verified experimentally. 

This work is supported by CINECA (PEC grant 2008) 
and MIUR (COFIN 2007). We acknowledge useful dis- 
cussions with A. Morgante, F. Becca and F. Mauri. 



[7; 

[8 

[9 
[10 

[11 

[12 

[13 

[14; 
[15 

[16 

[17 



[18 
[19 

[20 



[21 

[22 
[23 

[24 

[25 
[26 



* sorella@sissa.it 

] K. S. Novoselov et al., Science 306, 666 (2004). 

] L. Pauling, The nature of the chemical bond, 3rd ed. 

(Cor- nell University Press, Ithaca, NY, 1960), p. 204. 
.] A. K. Gcim, and K. S. Novoselov, Nature Mat. 6, 183 

(2007). 

] E. Rotenberg et al, Nature Mat. 7, 258 (2008); S. Y. 

Zhou et al., Natur e Mat. 7, 260 (200 8). 
.] Y.-M. Lu, Y. Ran. larXiv: 1005.4229k ! [cond-mat.str-el]. 
.] A. H. Castro Neto, Physics 2, 30 (2009); J. E. Drut and 
T. A. Lhde, Phys. Rev. B 79, 165425 (2009), and Phys. 
Rev. Lett. 102, 026802 (2009); W. Armour, S. Hands, 
and C. Strouthos, Phys. Rev. B 81, 125105, (2010). 
'] S. Pathak, V. B. Shenoy, and G. Baskaran, Phys. Rev. B 

81, 085431 (2010). 
.] A. H. Castro Neto et al., Rev. Mod. Phys. 81, 109 (2009). 
'] Z. Y. Meng et a l., Nature 464, 847 (2010). 
1] S. Sorella et al., |http://qe-f orge.org/projects/t urborvb/| 
] L. Spanu, S. Sorella, and G. Galli, Phys. Rev. Lett. 103, 
196401 (2009). 

] B. K. Clark, D. A. Abanin, and S. L. Sondhi, 

arXiv:1010.30TTV l [cond-mat.str-el]. 
■] M. Casula and S. Sorella, J. Chem. Phys. 119, 6500 
(2003). 

] M. Marchi et al., J. Chem. Phys. 131, 154116 (2009). 
] P. Fazekas and P.W. Anderson, Philos. Mag. 30, 423 

(1974); P.W. Anderson, Science 235, 1196 (1987). 
.] M. Burkatzki, C. Filippi, and M. Dolg, J. Chem. Phys. 

126, 234105 (2007). 
'] In particular our parametrization allows us to represent 
also slowly decaying pairing functions, as for instance 
the one obtained within the Hartree Fock theory, namely 
f(\r — r'|) oc \r — r'\~ 2 . 
] M. Capello et al., Phys. Rev. B 77, 144517 (2008). 
'] S. Azadi, C. Cavazzoni, and S. Sorella, Phys. Rev. B 82, 
125112 (2010). 

1] M. Casula, C. Filippi, and S. Sorella, Phys. Rev. Lett. 
95, 100201(2005); M. Casula et al., J. Chem. Phys. 132, 
154113 (2010). 

] J. Toulouse, and C. J. Umrigar, J. Chem. Phys. 128, 
174101 (2008). 

] C. J. Umrigar et al., Phys. Rev. Lett. 98, 110201 (2007). 
■] L. Bytautas, and K. Ruedenberg, J. Chem. Phys. 122, 
154110 (2005). 

] S. Parthiban, and J. M. L. Martin, J. Chem. Phys. 115, 
2051 (2001). 

.] B. Bernu et al., Phys. Rev. B 50, 10048 (1994). 
■] A. M. Black-Schaffer and S. Doniach, Phys. Rev. B 75, 
134512 (2007). 



